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TECHNICAL MEMORANDUM X-64601 

ANALYSIS OF A CONTROL SYSTEM 
CONTAINING TWO NONLINEARITIES 

SUMMARY 


A control system containing two nonlinear elements may be analyzed 
in the parameter plane to determine the existence, stability, and values of 
amplitude and frequency of limit cycles. By using describing functions to 
represent the nonlinearities, the system characteristic equation may be ob- 
tained. Two adjustable parameters are selected, each containing one of the 
describing functions and possibly several control system gains. A correla- 
tion between these parameters and the roots of the characteristic equation is 
determined by mapping stability contours from the complex s-plane onto the 
chosen parameter plane. A relationship between the inputs to the two nonlinear^ 
ities is determined next, and a locus representing the variation of the describ- 
ing functions also is plotted on the parameter plane . If this locus and the 
stability contour associated with a pair of pure imaginary roots intersect, the 
existence and characteristics of a limit cycle are indicated. Particular empha- 
sis is placed on organizing the equations so that they may be readily solved with 
a digital computer. 


INTRODUCTION 


A technique has been developed by Siljak [i ] to analyze systems con- 
taining two nonlinear elements if they have a common input or if their inputs 
are related by a linear differential equation. Viswanadham and Deekshatulu 
[2] and Gelb and Vander Velde [3] present techniques for analyzing a system 
whose inputs to the nonlinear elements are related by a nonlinear differential 
equation. A technique similar to that developed by Siljak is presented here to 
enable one to analyze this latter system. 

By using describing functions to represent the nonlinearities, the system 
characteristic equation maybe obtained. Two adjustable parameters are se- 
lected, each containing one of the describing functions and possibly several Con- 
trol system gains. A correlation between these parameters and the roots of 
the characteristic equation is determined by mapping stability contours from 
the complex s-plan onto the chosen parameter plane. 


A relationship between the inputs to the two nonlinearities is determined 
next, and a locus representing the variation of the describing functions also is 
plotted on the parameter plane. If this locus and the stability contour associated 
with a pair of pure imaginary roots intersect, the existence and characteristics 
of a limit cycle are indicated. From this point of intersection, the frequency 
and magnitude of the indicated limit cycle may be determined as a function of 
the characteristics of the nonlinearities and of the adjustable control system 
gains. The stability of the limit cycle is investigated by determining if all 
characteristic equation roots, other than the pair of pure imaginary roots, 
have negative real parts. This condition is indicated readily on the parameter 
plane. The behavior of the limit cycle when a small perturbation is applied to 
its amplitude also is apparent on the parameter plane. 

To demonstrate the application of this technique, a typical attitude con- 
trol system for a large space vehicle is analyzed. The results of the analysis 
are confirmed by analog simulation. 

ANALYSIS 

Hie system under consideration is portrayed in Figure i . It is assumed 
that the inputs Xj and X 2 to the nonlinearities are sinusoidal and that the applica- 
bility conditions permitting the use of the describing functions N t (A t , co) and 
N 2 (A 2 , ci) in lieu of the actual nonlinearities n^ and n 2 are satisfied, where 

x A = Aj sin cot 

and 

x 2 = A 2 sin ( out - ip) 

It also is assumed that only the first harmonic of the signal x^ passes through 
the first nonlinearity n^. Finally, it is assumed that the system parameters 
are time invariant. 

The describing functions are used to linearize the system characteristic 
equations; i.e. 

1 +N t Gj (s) N 2 G 2 (s) = 0, 
which leads to 

A < s ) = Z / R S k -0 . (2) 

k= 0 
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Figure 1. Control system containing two nonlinearities. 


where 

4 =b k p ‘ + °k p ’ +d k • < 3) 

Parameters and P 2 are defined as 

Pi = l t N t and P 2 = li N 2 , (4) 

where b^, c^, d^, lj, and 1 2 are constants or system parameters of the linear 
portion of the system . 

A correlation between parameters Pj and P 2 and the roots of the charac- 
teristic equation is determined by mapping certain stability contours from the 
complex s-plane onto the Pj-P 2 parameter plane. The stability contour chosen 
for this portion of the analysis is the imaginary axis of the s-plane. 

s = ift , (5) 

which will be referred to as the t, = 0 contour. To facilitate computation and 

to provide flexibility for use in later analysis, s of equation (2) maybe re- 
placed by 

sk "\ +iY k • (6) 

where and may be obtained by using the recurrence formulas [1 ]: 
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and 


*k + i + *“n t ^ + “n , Vl-° 


Y +2w tY, + i*> 2 Y, , =0. 
k+i n k n k-i 


(7) 


Equation (2) may be separated into its real and imaginary parts and regrouped: 

n 


Re{A (s)} = Sf X = Bi Pj + Ci P 2 + D t = 0 
k-0 k 


(8. a) 


and 


n 

Im{A(s)} - 2 4 \ =B 2 p l +C 2 +^2 =0 . 


k-0 


(8.b) 


where 


B * - 2 Vk 

k= 0 

B; = 2 b Y 
k=o k 


n n 

Cl = 2 °|X k • D, = 2 \X,C • 


k-0 

n 


k-0 


n n 

c 2 = 2 c k Y k . Di = 2 y k ■ 


k-0 


k= 0 


Equations (8. a) and (8.b) maybe solved for Pj and P 2 which may be plotted on 
the P 1 -P 2 parameter plane as functions of frequency ft. This curve is the 
= 0 stability contour associated with complex conjugate roots. The boundary 
separating stable real roots from unstable ones is determined by setting s = 0 
in equation (2) and solving for P 2 as a function of Pj. The resulting curve also 
is plotted on the P t -P 2 parameter plane. The stable region (corresponding to 
the left half of the s-plane) in the parameter plane is determined by observing 
the sign of the Jacobian (J) [1]. 


In determining regions of stability and instability on the parameter plane, 
one must be alert for the possible existence of singular cases where the 
Jacobian J = 0 2 . This may happen for certain combinations of and N 2 — a 
point not mentioned by Viswanadham and Deekshatulu [2]. (Analysis of the 
Jacobian associated with the example presented in Viswanadham's and 
Deekshatulu's paper [2] shows the existence of several singular points for cer- 
tain combinations of values of Nj, N 2 , and u>.) In such cases, the mapping of 
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(, 


certain points from the complex s -plane onto the parameter plane may result 
in loci of points or in distinct points. Hence, additional stability boundaries 
may occur on the parameter plane. One such singular/case corresponds to the 
mapping of the origin of the s-plane onto the parameter plane as the real root 
stability boundary. 

To obtain some insight into the choice of values for adjustable gains of 
the linear portion of the system, two such gains (or combinations of gains) may 
be chosen as parameters ( o \ and o^) • Considering linear operation of the 
system, the characteristic equation (2) may be written in terms of and 
by setting Nj and N 2 equal to unity; i.e. , 

4 =b k “ 1+C k c%+d k • <9) 

Stability contours associated with chosen values of the damping ratio £, may be 
plotted as functions of on the Oj-Qfc parameter plane by mapping the line 

s = - £C0 + id) ^( i - £*) (10) 

n n 

from the s-plane onto the a^-a^ parameter plane. Any roots lying on this line 
will possess an associated damping ratio specified by £. If equation (9) is 
substituted for equation (3) and P t and P 2 are replaced by and afc in equa- 
tions (8. a) and (8.b) , they may be solved by and for various values of 

£ . Then, values of a.y and a» may be chosen for desired values of £ and w . 

n 

The relationship between and x 2 may be determined, assuming that 
Gi(s) is a low-pass function and that only the first harmonic passes through n t . 
Since N* (Aj, cu) is known, A 2 may be calculated as a function of Aj and co : 

A 2 = |G t (io>) |N| (A,, w) A t . (11) 

For chosen values of A l and oi, N 2 may be found using equation (11) . Thus, N 2 
may be plotted versus Nj as a function of Aj and co. Using equation (4) , one may 
plot this nonlinear locus on the Pi-P 2 parameter plane, as a function of A lt for 
various values of cu. If the locus intersects the £ = 0 curve, a limit cycle is 
indicated if the frequencies of the locus (w) and the £ = 0 curve ( ft) match; 
this is the frequency of oscillation of the limit cycle. The amplitude of the 
limit cycle is determined from the value of Aj corresponding to the intersection. 
For a limit c^cle"to“be‘ stable ,'it _ is necessary but not sufficient -that all-roots- - ~ 
(other than the pair of complex roots lying on the imaginary axis of the s-plane) 
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lie to the left of the imaginary axis. Hence, the intersection of the locus and 
the f = 0 curve must lie on a boundary of the stable region ( the entire f = 0 
curve is not necessarily a stability boundary — a point omitted by 
Siljak [1]) • If the locus crosses from the unstable region as Aj increases, a 
stable limit cycle is indicated (and vice versa if it crosses from the stable 
region into the unstable region) . If no intersection occurs and the locus remains 
in the stable region for all values of A lf linear stable operation is indicated. If 
the locus remains in the unstable region, unstable operation without a limit 
cycle is indicated. 

Example. Consider a typical pitch -plane attitude control system for a 
space vehicle (Fig. 2) . The rigid body dynamics are represented by the trans- 
fer function 

G KB (S ) S ® (s)//3 l ( s) = c/s , (12) 

where 

c = (BL/Bp- )(3 r 1 /I . 

Li Lg 


The distance between the center of mass and the point of application of the 

thrust is denoted by 1 , the lateral component of thrust by L, and the moment 

S 

of inertia about the pitch axis by I. 


The system characteristic equation may be written in the form of 
equation (2) : 

3 

A(s) — t 3 + s 2 +a t cN N s + a 0 cN = 2 fi ~ 0 > 

s s u s 

k= 0 


(13) 


where and represent describing functions associated with the saturation 

characteristic and the dead-zone characteristic, respectively. Convenient 
parameters are 

Po = a 0 N and Pj s a t N N , (14) 

s s u 


where 





(15) 
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CONTROL SERVOACTUATOR 



-.THRUST VECTOR ATTiTMOC 

Figure 2. Space vehicle attitude control model. 

In the manner indicated by equations (8. a) and (8.b) , the £ = 0 contour may be 
determined: 

P 0 = nVc (16) 

and 

Pi = rfiVc = tP 0 , (17) 

where 

J =c 2 fi > 0, Vfi >0 . (18) 

In a design problem it is desirable to separate the describing functions N and 

s 

Np rather than have them appear as a product in the parameter Pi. This is 
achieved by defining P t : 

Pi =ai Np/a 0 . (19) 

Examination of the Jacobian reveals that it will become identically equal to 
zero only if o> equals zero ( since c is non-zero) . This corresponds to the sin- 
gular case associated with the real root stability boundary which is 

P 0 = 0. (20) 
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Now the stability boundaries may be plotted, using equations ( 16) , ( 19) , and 
(20) (Fig. 3) . The number of stable roots of the characteristic equation is 
indicated within square brackets for each region. 

To determine suitable values for the adjustable control gains, a 0 and a t , 
linear operation is assumed and contours corresponding to various values of £ 
maybe plotted as functions of ok. A typical contour (for £ = 0.5) also is 

shown in Figure 3; values of o> and fi are shown in parentheses. Numerical 

n 

values must be assigned to the constants r, a 0 , a*, k , k , S, D, and c. 

S JJ 

Constants r, k^, and D are fixed by the design of the actuator. Values repre- 
sentative of contemporary hardware are r = 0.2s, k g =1, and S = 1 deg. A 
typical rate gyroscope has values of k^ = 1 and D = 0. 1 deg/s. Parameter c is 

a function of flight time and, for a typical launch vehicle, varies slowly between 
0.5 and 2.5. The value of 0.7 is chosen from the example. Values of a 0 and aj 
are chosen from Figure 3, assuming linear operation and £“0.7 and u » 0.6 
rad/s; a 0 =0.5 and a t = 1. 

Assuming the following forms for x* and x 2 : 

Xi = /3 = A sin fit 
1 s 

and (2i) 

Xg = 3 = A^ sin ( fit + ipi) 


and utilizing the assumptions of the foregoing section, the relationship, between 
Ap and A g is 


A^ =cN (A , w) A /o> . 
D s s s 


( 22 ) 


Describing function relationships for N (k , A , S) and N (k_^, A^, D) are 

s s s D D D 

well known and readily available. From these relationships and equations (16), 

(17) , (20) , and (22) , loci representing the variation of N and N may be 

s u 

plotted on the Pn-Pi parameter plane as a function of A for various values of co. 

s 

These loci are shown as dashed lines in Figure 3. Sets of two figures within 

parentheses indicate corresponding values of A /S and A /D. Since most loci 

s u 

shown rise vertically at P 0 =0.5, the values shown for A g /S and A^/D 
correspond to the locus associated with u> = 0.592. It is seen that a limit cycle 
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Figure 3. P 0 -Pi parameter plane 
(a 0 =0.5, a t = 1, r = 0.2, and 
c = 0.7). 


cannot be avoided if otherwise nonlinear 
operation is to be precluded . It also 
is observed that only one limit cycle is 
possible, and it will be stable. The 
characteristics of the limit cycle for 
this example are read from the inter- 
section of the f = 0 contour and the 
nonlinear locus as: fi =0.592; 

P 0 =0.5; Pi =0.2; A =0.106, S = 1; 

s 

and Aq = 0.126, D =0.1. The predicted 

numerical results were confirmed by 
simulating the system on an analog 
computer, thereby increasing confi- 
dence in the technique and the assump- 
tions. Figures 4(a) and 4(b) show 
these results for several initial condi 4 
tions on 0. Also of interest (although 
not analyzed) is the case where an 
input signal is applied. Simulation 
results, shown in Figure 5, indicate 
that limit cycle operation will not 
occur when the system is forced by a 
ramp input, but as soon as the input 
levels off to a step input, sustained 
oscillations again occur. 


In the actual design of a control system, the values chosen for t, a 0 , 
a l, k , k , S, and D would be varied (within constraints imposed by actual 

S A/ 

hardware implementation) to achieve other limit cycle characteristics so that 
the least undesirable characteristics could be obtained. Also, the problem 
would be analyzed for the entire range of values that c assumes during predic- 
ted flight . 

CONCLUSIONS 

The parameter plane technique of stability analysis can be applied to a 

model of a system containing two nonlinearities whose inputs are related by a 
nonlinear differential equation. The constraints imposed on the system are 
prescribed, and limit cycle operation and characteristics may be predicted. 

In cases where limit cycle operation can be avoided, relative stability 
characteristics can be established by choosing the distribution of the roots of 
the characteristic equation through the variation of two adjustable parameters. 
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Figure 4(a) . Computer simulation: 
0 (0) =0.4 deg (time scale: 2.5 


Figure 4(b) . Computer simulation: 
0 (0) =0.1 deg (time scale: 2.5 
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Figure 5. Computer simulation: forced solution (time scale: 2.5 s/div) 
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